/**
 * This file is part of gapro library.
 * 
 * gapro library is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
 * later version.
 * 
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
 * details.
 * 
 * You should have received a copy of the GNU Lesser General Public License along with this program. If not, see
 * <http://www.gnu.org/licenses/>
 * 
 * Copyright 2012 Vinay Penmatsa. All rights reserved.
 * 
 */
package org.stupidwhiz.gapro.core;

/**
 * Matrix operations with Fraction types
 * 
 * @author Vinay Penmatsa
 * @since 1.0
 */
public class MatrixFraction {

    private static final String IDENTITY_AS_STRING = "<1,0,0> <0,1,0> <0,0,1>";
    Fraction[][] matrix;

    public MatrixFraction(Fraction[][] matrix) {
        this.matrix = matrix;
    }

    /**
     * Returns the number of rows
     * 
     * @return row count
     */
    public int rowCount() {
        return matrix[0].length;
    }

    /**
     * Return the number of columns
     * 
     * @return column count
     */
    public int columnCount() {
        return matrix.length;
    }

    /**
     * Returns a new matrix that is the transpose of this matrix.
     * 
     * @return the tranpose of this matrix
     */
    public MatrixFraction transpose() {
        int m = columnCount();
        int n = rowCount();
        Fraction[][] transpose = new Fraction[n][m];

        for(int i = 0; i < n; i++) {
            for(int j = 0; j < m; j++) {
                transpose[i][j] = matrix[j][i];
            }
        }

        return new MatrixFraction(transpose);
    }

    /**
     * Returns a new matrix that is the inverse of this matrix
     * 
     * @return the inverse of the matrix
     */
    public MatrixFraction inverse() throws Exception {
        int m = columnCount();
        int n = rowCount();

        if(m != n) {
            throw new Exception("Matrix is not square");
        }

        MatrixFraction im = identity(m);
        MatrixFraction augmented = augment(this, im);
        reduce(augmented);

        MatrixFraction id = subMatrix(augmented, 0, 0, m, m);
        if(!isIdentity(id))
            throw new Exception("Matrix not invertible");

        return subMatrix(augmented, 0, m, n, m);
    }

    private static boolean isIdentity(MatrixFraction mf) {
        return IDENTITY_AS_STRING.equals(mf.toString());
    }

    /**
     * Returns a sub matrix starting at row, col and of dimensions rowLength and colLength
     * 
     * @param row
     * @param col
     * @param rowLength
     * @param colLength
     * @return
     */
    public static MatrixFraction subMatrix(MatrixFraction mat, int row, int col, int rowLength, int colLength) {
        if(mat == null || row < 0 || col < 0 || rowLength < 0 || colLength < 0)
            throw new IllegalArgumentException();

        int m = mat.columnCount();
        int n = mat.rowCount();

        if(row + rowLength > n || col + colLength > m) {
            throw new IllegalArgumentException("Requested dimension is greater that size of matrix");
        }

        Fraction[][] f = new Fraction[colLength][rowLength];

        for(int j = 0; j < colLength; j++) {
            for(int i = 0; i < rowLength; i++) {
                f[j][i] = mat.matrix[col + j][row + i];
            }
        }
        MatrixFraction mf = new MatrixFraction(f);
        return mf;
    }

    private void reduce(MatrixFraction mx) {
        int row = 0; // row index
        // loop columns 0 to colLength-1

        for(int col = 0; col < mx.columnCount(); col++) {
            int k = findMaxInRowFromIndex(mx.matrix[row], row);
            if(k == -1)
                continue;
            if(k != row) {
                exchangeRows(mx, k, row);
            }
            Fraction m = mx.matrix[col][row];
            multiplyRow(mx, row, m);
            reduceRow(mx, row, col);
            row++;
            if(row == mx.rowCount())
                break;
        }

    }

    private int findMaxInRowFromIndex(Fraction[] row, int index) {
        Fraction max = row[index];
        int k = index;
        for(int i = index; i < row.length; i++) {
            if(row[i].compareToAbs(max) > 0) {
                max = row[i];
                k = i;
            }
        }
        if(max.getNumerator() == 0)
            k = -1;
        return k;
    }

    private void exchangeRows(MatrixFraction mx, int k, int row) {
        for(int j = 0; j < mx.columnCount(); j++) { // col
            Fraction temp = mx.matrix[j][k];
            mx.matrix[j][k] = mx.matrix[j][row];
            mx.matrix[j][row] = temp;
        }
    }

    private void multiplyRow(MatrixFraction mx, int row, Fraction m) {
        Fraction invM = new Fraction(m.getDenominator(), m.getNumerator());
        for(int i = 0; i < mx.columnCount(); i++) {
            mx.matrix[i][row].multiply(invM);
        }
    }

    private void reduceRow(MatrixFraction mx, int row, int col) {
        for(int r = 0; r < mx.rowCount(); r++) {
            if(r == row)
                continue;
            Fraction f1 = new Fraction(mx.matrix[col][r].getNumerator(), mx.matrix[col][r].getDenominator());
            for(int c = 0; c < mx.columnCount(); c++) {
                Fraction f = new Fraction(f1.getNumerator(), f1.getDenominator());
                f.negate();
                f.multiply(mx.matrix[c][row]);
                mx.matrix[c][r].add(f);
            }
        }
    }

    public static MatrixFraction augment(MatrixFraction m, MatrixFraction n) throws Exception {
        int ml = m.rowCount();
        int nl = n.rowCount();

        if(ml != nl) {
            throw new Exception("Row count of m, n must be equal");
        }

        int colCount = m.columnCount() + n.columnCount();
        Fraction[][] augmented = new Fraction[colCount][ml];

        for(int i = 0; i < m.columnCount(); i++) {
            System.arraycopy(m.matrix[i], 0, augmented[i], 0, ml);
        }

        for(int i = 0; i < n.columnCount(); i++) {
            System.arraycopy(n.matrix[i], 0, augmented[i + m.columnCount()], 0, nl);
        }

        return new MatrixFraction(augmented);
    }

    /**
     * Creates a n x n identity matrix. Elements are initialized to 0 and diagonal elements to 1
     * 
     * @param n
     * @return identity matrix
     */
    public static MatrixFraction identity(int n) {
        Fraction[][] m = new Fraction[n][n];
        int key = 0;
        for(int d = 0; d < n; d++) {
            for(int f = 0; f < n; f++)
                m[d][f] = new Fraction(0);
            m[d][key] = new Fraction(1);
            key++;
        }
        return new MatrixFraction(m);
    }

    /**
     * Returns the string representation of the matrix as column vectors
     * 
     */
    public String toString() {
        StringBuilder sb = new StringBuilder();
        for(int j = 0, m = columnCount(); j < m; j++) {
            sb.append("<");
            for(int i = 0; i < rowCount(); i++) {
                sb.append(matrix[j][i] + (i == rowCount() - 1 ? "" : ","));
            }
            sb.append(">");
            if(j != m - 1)
                sb.append(" ");
        }
        return sb.toString();
    }
}
